Diverse electrical responses in a network of fractional-order conductance-based excitable Morris-Lecar systems

The diverse excitabilities of cells often produce various spiking-bursting oscillations that are found in the neural system. We establish the ability of a fractional-order excitable neuron model with Caputo’s fractional derivative to analyze the effects of its dynamics on the spike train features observed in our results. The significance of this generalization relies on a theoretical framework of the model in which memory and hereditary properties are considered. Employing the fractional exponent, we first provide information about the variations in electrical activities. We deal with the 2D class I and class II excitable Morris-Lecar (M-L) neuron models that show the alternation of spiking and bursting features including MMOs & MMBOs of an uncoupled fractional-order neuron. We then extend the study with the 3D slow-fast M-L model in the fractional domain. The considered approach establishes a way to describe various characteristics similarities between fractional-order and classical integer-order dynamics. Using the stability and bifurcation analysis, we discuss different parameter spaces where the quiescent state emerges in uncoupled neurons. We show the characteristics consistent with the analytical results. Next, the Erdös-Rényi network of desynchronized mixed neurons (oscillatory and excitable) is constructed that is coupled through membrane voltage. It can generate complex firing activities where quiescent neurons start to fire. Furthermore, we have shown that increasing coupling can create cluster synchronization, and eventually it can enable the network to fire in unison. Based on cluster synchronization, we develop a reduced-order model which can capture the activities of the entire network. Our results reveal that the effect of fractional-order depends on the synaptic connectivity and the memory trace of the system. Additionally, the dynamics captures spike frequency adaptation and spike latency that occur over multiple timescales as the effects of fractional derivative, which has been observed in neural computation.

www.nature.com/scientificreports/ generate (depending on fractional order exponents) a wide range of firing phenomena or multiple timescale dynamics [10][11][12][13][14][15][16][17][18][19][20][21] . As such, a deeper understanding has been reached in different areas of biophysical processes 11,[22][23][24][25][26][27][28] , showing more realistic dynamical features [29][30][31][32][33] . Fractional-order derivatives provide a mathematical framework in which memory dependent properties are considered. Earlier, a geometrical 7,34 interpretation for the fractionalorder derivative was introduced, which suggests that inhomogeneity of the time scale exists in the system. It may have an impact on the delays of signals or history dependent activities in comparison to the temporal order dynamics. In the fractional-order domain, the present state of the system is influenced by the previous states. History dependent spiking features are important, as the neuronal activities develop over time and continuously integrate the previous information 9, 35 . Fractional-order neuronal models have been applied to study different firing responses 13,18,19,36,37 , firing rates and spike frequency adaptation both theoretically and experimentally. For some fractional-orders less than one, initially, the voltage increases faster, however, it reaches the steady state condition, i.e., quiescent state after longer time duration. Neurons can show adaptation in the fractional domain when we scale the input stimulus. The single spikes and various bursting maintain different information. The adaptation depends on the fractional exponent and the mean firing rates can change and adapt to the variations in the stimulus. The spike frequency adaptation follows power law dynamics in the fractional-order domain 10,18,38 .
The primary goal of the paper is to provide a brief description in understanding the effects of fractionalorder derivative on the electrical activities of single M-L spiking neurons with class I & class II excitabilities and the slow-fast M-L neurons 1,2,39 with its network architecture. One approach to study such firing characteristics and adaptation is to consider a conductance-based model that explores the intrinsic dynamics underlying the fractional-order derivatives. Previous works have investigated the various spiking responses depending on various parameter regimes, however, FOD can itself explore diverse firing responses 13,16 . The M-L models are taken into account for their diverse responses ranging from spiking to bursting. We consider different regimes in the parameter space of the M-L model: tonic spiking and fast spiking. Further, the model is extended to its 3D counterpart, where the applied current, I is not constant, but rather varies with time. We consider the fractionalorder as the predominant parameter in the system and when it changes slowly, the spike transitions occur and we observe mixed-mode oscillations (MMOs) and mixed-mode bursting oscillations (MMBOs). It is one of the most interesting neuronal oscillations that emerge from the electrical activities 40,41 . MMOs are used to describe the alternating trajectories between small and large amplitude oscillations (SAOs and LAOs) 42,43 . These make the system fascinating and the output provides interesting and potential applications in a dynamical system. The emergence of MMBOs creates a spike adding mechanism. Earlier, it was observed that the MMOs reviewed the dynamical and neuronal behavior of locomotion or breathing 44,45 . It was observed in calcium signaling and electrocardiac systems 46,47 . Krupa et al. 48 examined the mechanism of MMOs oscillations in a two-compartmental model of dopaminergic neurons in the mammalian brain stem. We also investigate the impact of electrical coupling on a mixed population where the neurons are either quiescent or oscillatory. Here, the neurons are assumed to be connected through the links of the Erdös-Rényi network. The coupling induces complex firing activities such as periodic bursting or spike frequency adaptation for all the nodes in the network. Based on the observed synchronization phenomena, a reduced-order model is developed which can produce the activities of the entire network.
In our work, we find consistent differences in the characteristics of the neuronal functional behavior using the fractional exponent. The fractional-order voltage dynamics can significantly change the spiking features of different single neuron models [12][13][14][15]17,18,36,37 . Realistic features can build the model more sensitive to neuronal dynamics, particularly in the potential collective behavior of the network, where past dynamical behavior might influence the present states.

Formulation of the excitable model dynamics and some preliminaries
In this section, we describe the fractional-order excitable conductance-based model and review the existence of various characteristics observed in cortical areas 2 . We establish a particular parameter regime that supports the firing features with the variations of fractional exponent. To generate diverse spikes using fractional-order dynamics, we study the 2D and 3D M-L models with particular parameters and channel dynamics. Here, we choose the two models to separate the effects of fractional derivatives on the dynamical behavior of the model. Morris and Lecar 1,2 proposed a simple mathematical model to describe the oscillations in the barnacle giant muscle fiber consisting of the membrane voltage equation with instantaneous activation of calcium current and an additional recovery equation describing slow activation of potassium current. The 2D M-L model is described in a commensurate fractional-order domain as follows The biophysically motivated excitable model involves a voltage-gated Ca 2+ current, delayed rectifier K + current and a leak current respectively. u measures the membrane voltage dynamics and v is the activation variable of K + ion channels. The parameters g Ca , g K and g L indicate the maximum conductance functions to Ca 2+ , K + and leak currents respectively. V Ca , V K and V L are the reversal potentials to different ionic current functions. C measures the membrane capacitance. φ represents the temperature scaling factor for K + channel opening. The parameters V 1 , V 2 , V 3 and V 4 have fixed positive values. I indicates the applied stimulus. We would like to account the effects of various injected current stimulus on the fractional-order system with the fractional exponent, α ( 0 < α ≤ 1).
(1) www.nature.com/scientificreports/ Slow-fast dynamical phenomenon. First, we assume the neuron is at the onset of firing and it generates spike generation as the control parameter moves slowly. The slow-fast dynamics can be mathematically modeled as 1 where ẋ(t) = f (x, z) (fast spiking) and ż(t) = δg(x, z) (slow modulation). x ∈ R m represents the fast variables and z ∈ R n the slow variables with 0 < δ << 1 measuring the timescale separation parameter.
The following system of ODEs represents the slow-fast 3D M-L model where (u, v) denote the fast subsystem and w slow variable. The fractional-order modified 3D M-L model 1,13 is presented as The system has the following characteristics. The system variable w is the external injected current which follows the power law dynamics in the fractional system and characterizes the memory effect of the membrane potential 1,13,16 . The parameters V 1 , V 2 , V 3 and V 4 are suitably selected for the hyperbolic functions in order to explain that they can reach their equilibrium points instantaneously. The parameter value µ is less than 1 i.e., 0 < µ < 1 which measures the ratio of time scale between oscillations and modulation. Lundstrom et al. 10 studied that pyramidal neurons can act as fractional differentiators of the stimulus amplitude envelope for this type of input. FOD can generalize the derivative operator such that, to obtain the first order derivative of a function, differentiate twice taking the fractional-order derivative of order α = 1/2 and it results in the first derivative 4,6 . It filters the response with a decaying time constant that depends on α.
The parameter sets for all the simulation results are considered as (for Eq. 1) 1,39 Set I: C = 20, g Ca = 4, .067 (for class I excitable membrane model) with varying I, Set I: I = 40 , Set II: I = 45 and for the class II membrane model, the parameters are the same described above except Set III: g Ca = 4.4, V 3 = 2, V 4 = 30, φ = 0.04 , and I = 100 . In order to study the system dynamics, we first analyze the equilibrium states and then bifurcations. Next, we use the following sets of parameters to deal with system (3) and its modified versions by considering I(w) = 0.08 − 0.03w using C = 1 for the parameter sets I, II and III respectively.

Preliminaries to systems of fractional-order differential equations
To study the fractional-order M-L model, we consider the familiar definition of the fractional derivative i.e., the Caputo fractional-order derivative 6,34 . The commensurate fractional-order model with fractional exponent α ∈ (0, 1) can be described as for 2D and 3D cases, respectively. The Caputo fractional differential operator is defined as where the Gamma function is given by Ŵ(z) = ∞ 0 e −s s z−1 ds . The limits of the integration (i.e., from 0 to t) show that, in contrast with the classical integer-order derivative, the fractional-order derivative depends on the whole previous history of the function. Hence, due to the non-locality of the Caputo differential operator, a fractionalorder mathematical model is able to reflect memory properties of the system variables. It is important to note that for α = 1 , the Caputo derivative converges to the first-order integer derivative. An additional advantage of Caputo-type fractional -order derivative over other types of fractional differential operators is that the derivative of a constant is zero. It is efficient to integrate all the previous activities of the function weighted by a function that follows power-law dynamics 9,14,15 .

Remark 3.1
In the investigation of the local stability properties of an equilibrium of a dynamical system, the classical Hartman-Grobman linearization theorem plays a fundamental role: it states that the local behavior of www.nature.com/scientificreports/ a dynamical system in a neighborhood of a hyperbolic equilibrium is qualitatively equivalent to the behavior of its linearization near the equilibrium. It is important to remark that a fractional-order counterpart of this linearization theorem has been obtained in 30 . If X ⋆ is an equilibrium of system system (4), i.e. f (X ⋆ ) = 0 , the corresponding linearized system at X ⋆ is: where J f (X ⋆ ) is the Jacobian matrix of the function f computed at X ⋆ . Therefore, the equilibrium X ⋆ of the nonlinear system (4) is asymptotically stable if and only if the trivial solution of the linearized system (6) is asymptotically stable [49][50][51] . Furthermore, based on the well-known Matignon's theorem 49 , the linearized fractional-order system (6) is asymptotically stable if and only if | arg( )| > απ 2 , for any eigenvalue of the Jacobian matrix J f (X ⋆ ).

Definition 3.1
If some eigenvalues of the Jacobian matrix J f (X ⋆ ) satisfy | arg( )| > απ 2 and some other eigenvalues satisfy | arg( )| < απ 2 , then the equilibrium X ⋆ is a called a saddle point 52,53 .

Remark 3.2
In a 3D nonlinear fractional-order system, an equilibrium X ⋆ is called a saddle of index one if one of the eigenvalues of the Jacobian matrix J f (X ⋆ ) is unstable (i.e. | arg( 1 )| < απ 2 ) and other two eigenvalues are stable | arg( 2,3 )| > απ 2 . On the other hand, if two eigenvalues associated to the equilibrium X ⋆ are unstable, while only one eigenvalue is stable, the saddle point X ⋆ is called saddle of index two 53 .
We numerically simulated the model (Eqs. 1 and 3) using the L1 scheme 4,12,18 and approximated the fractionalorder derivative as described in Appendix.

Qualitative analysis and theoretical considerations
Analysis of the 2D system. System (1) is a particular case of the 2D fractional-order conductance-based excitable model: where u and v are the membrane voltage and the gating variable of the neuron, I is an applied current, Ĩ (u, v) represents the ionic current, ℓ(v) is the rate constant for opening ionic channels and v ∞ (v) represents the fraction of open ionic channels at steady state, respectively.
In particular, we have from model (1): and The equilibrium points of system (7) are the solutions of the algebraic system: which is equivalent to The function I ∞ (u) satisfies the following basic properties: has exactly two real roots u max < u min .
Denoting by I max = I ∞ (u max ) and I min = I ∞ (u min ) the maximum and minimum values of I ∞ , respectively, the function I ∞ is increasing on the intervals (−∞, u max ] and [u min , ∞) and decreasing on the interval (u max , u min ). Hence, depending on the external input I, there are exactly three branches of equilibrium points, denoted by Consequently, one of the following situations may hold: • If I < I min or if I > I max , then system (1) has a unique equilibrium point.
. www.nature.com/scientificreports/ • If I = I min or if I = I max , then system (1) has two equilibrium points.
• If I ∈ (I min , I max ) , then system (1) has three equilibrium points.
The Jacobian matrix associated to the system (1) at an arbitrary equilibrium state In this case, the necessary and sufficient conditions for the asymptotic stability of an equilibrium point (u ⋆ , v ⋆ ) reduce to the following inequalities 54 : where We first remark that the second branch of equilibrium points is completely unstable. Indeed, any equilibrium point (u 2 (I), v ∞ (u 2 (I))) with I ∈ (I min , I max ) , satisfies I ′ ∞ (u 2 (I)) < 0 , and hence, δ(u 2 (I)) < 0 . In fact, the negative sign of the Jacobian's determinant guarantees that each equilibrium point of the second branch is a saddle point, no matter what fractional-order α is considered in system (7).
On the other hand, along the first and third branches of equilibrium points, it is easy to find that the determinant of Jacobian is positive. Hence, the stability of the equilibrium points particularly depends on the trace τ . Obviously, if τ (u ⋆ ) < 0 , the equilibrium point (u ⋆ , v ⋆ ) is asymptotically stable, irrespective of the fractional-order α considered in system (7). However, if τ (u ⋆ ) ≥ 0 , an equilibrium point (u ⋆ , v ⋆ ) of the first or the third branch is asymptotically stable, if and only if We will further assume that V k < u max < u min < 1 . We can easily evaluate: is an equilibrium point of the third branch such that u ⋆ > 1 , it follows that Ĩ u (u ⋆ , v ∞ (u ⋆ )) > 0 , and hence τ (u ⋆ ) < 0.
Furthermore, we can also express and hence, if (u ⋆ , v ⋆ ) is an equilibrium point of the first branch such that u ⋆ < V K , we deduce that Ĩ u (u ⋆ , v ∞ (u ⋆ )) > 0 , and similarly as above, we get τ (u ⋆ ) < 0. Based on the above calculation, we also remark that: for either u m = u max or u m = u min , and assuming that φ is small enough, it can be observed that the inequality τ (u m ) > 0 might hold. Therefore, the function τ (u) might have two roots u ′ ∈ (V K , u max ) and u ′′ ∈ (u min , 1) , respectively. Based on the numerical data, we can further assume that if they exist, these roots are unique in the aforementioned intervals.
In conclusion, the stability of equilibrium states may depend on the fractional-order α only in the following two cases: • the equilibrium point belongs to the first branch and u ⋆ ∈ (u ′ , u max ); • the equilibrium point belongs to the third branch and u ⋆ ∈ (u min , u ′′ ).
In this case, the critical value α ⋆ given by (9) corresponds to a Hopf-type bifurcation (i.e. the Jacobian matrix has a pair of complex conjugate eigenvalues such that | arg( )| = απ 2 ). In other words, the position of the Hopf bifurcation points in the (I, u)-plane, situated on the first and / or third branches, respectively, depending on the fractional-order α considered in system (7). Obviously, this will have a direct effect on the type of spiking and bursting behavior both in the 2D system (7) as well as in the 3D slow-fast system, as it will be unveiled in the next section.
As an example, first we show the bifurcation scenario of the classical 2D M-L model with Hopf bifurcation points by considering I as a bifurcation parameter (Fig. 1a). Next, we consider the effect of fractional order on the dynamics of system (1) and showed that how it stabilizes the system as we decrease the value of α (Fig. 1b). Then, Fig. 2 presents the phase portraits of the 2D M-L model with the parameters from Set II, for different www.nature.com/scientificreports/ values of the fractional-order α . In this case, there is only one unstable equilibrium point for the system, namely (u ⋆ , v ⋆ ) = (5.08955, 0.311245) (at the intersection of the nullclines), situated on the third branch. The critical value of the fractional-order corresponding to the Hopf bifurcation at the equilibrium point is α ⋆ = 0.787825 , computed by the formula (9). In the integer-order case, α = 1 , a large-amplitude limit cycle attractor is present, corresponding to spiking behavior. As the fractional-order α decreases, the large-amplitude attractive quasiperiodic limit cycle approaches the unstable equilibrium point and as α approaches the critical value α ⋆ for Hopf  www.nature.com/scientificreports/ bifurcation, a more complex quasi-periodic orbit emerges, involving smaller-amplitude oscillations around the equilibrium, as well as large-amplitude spikes. When α < α ⋆ , the equilibrium becomes asymptotically stable. In Fig. 3, we show corresponding time series to further verify the numerical results, noting that the computed critical values of the fractional-order are α ⋆ = 0.757245 for Set I and α ⋆ = 0.834537 for Set III, respectively.
Analysis of the 3D system. In line with the previously presented aspects, the 3D slow-fast fractionalorder model (3) can be written as: where Ĩ (u, v) is given by (8) and We will assume that I(w) and V 3 (w) are decreasing functions, and that V 0 + V K < 0 (according to the considered parameter sets). The unique equilibrium point of system (10) is and w ⋆ is the unique root of the strictly decreasing function w → I(w) −Ĩ(−V 0 ,ṽ(−V 0 , w)).
The Jacobian matrix at the equilibrium point (u ⋆ , v ⋆ , w ⋆ ) is and its charactersitic equation is: As c > 0 , it follows that the product of the eigenvalues of the Jacobian matrix is negative, and hence, one of the eigenvalues is a negative real number and the other two eigenvalues are either complex conjugated, or are real and have the same sign. We will further assume that at least one of the coefficients a or b is negative (based on the parameter sets under consideration), and hence, it is clear that the Routh-Hurwitz conditions are not satisfied for the characteristic polynomial. Hence, denoting by the discriminant of the characteristic polynomial, we distinguish two cases: • if � > 0 , the Jacobian matrix J has one negative and two positive eigenvalues, and consequently, the equilibrium point (u ⋆ , v ⋆ , w ⋆ ) is a saddle point of index two, for any fractional-order α (e.g. in the case of parameters from Set I and Set II); • if � < 0 , the Jacobian matrix J has one negative eigenvalue and two complex conjugate eigenvalues with positive real part (e.g. in the case of parameters from Set III). Consequently, there exists a critical value α ⋆ of the fractional-order such that the equilibrium point (u ⋆ , v ⋆ , w ⋆ ) is asymptotically stable for α < α ⋆ and unstable for α > α ⋆ . At α = α ⋆ , a Hopf-type bifurcation occurs in a neighborhood of the equilibrium point, resulting in the appearance of persistent oscillations. The critical value α ⋆ is found using the method presented in 55 ( α ⋆ = 0.62477 for Set III).

Analysis of diverse oscillatory responses
We start our discussion with the fractional-order class I and class II single M-L neurons and then extend it to the slow-fast dynamics. We simulated the spikes from the single model, and the membrane voltage dynamics depends on the voltage-gated conductances. The input stimulus is considered as I. We tuned the fractional-order exponent, α , with different parameter regimes: tonic spiking and fast spiking zone for class I and class II excitabilities. We next show the modulations of the electrical activities for a long time scale and the spike frequency adaptive effects. We considered two different suitable current stimuli, I = 40 and 45 for class I neuron and I = 100 for class II case. We choose these types of input stimuli as it shows the tonic spiking and fast spiking for the classical-order dynamics, however, when we change it in the fractional domain, the dynamical model produces variations in the firing features not explored earlier to the best of our knowledge. The bifurcation analysis is performed and the numerical results are supported by the stability analysis and there is a good agreement between the analytical and numerical findings. First, we consider the class I excitable M-L neuron with parameter sets I and II. The classical-order neuron shows tonic spiking when stimulated, when the input stimulus current is on ( I = 40 ), the neuron continues to exhibit a train of spikes, called tonic spiking. Then, as the fractional exponent decreases to α = 0.85 , it shows tonic spiking however, the interspike interval increases, i.e., firing frequency decreases. With further decrease of α = 0.83 and 0.80, it generates regular bursting and then regular bursting with low firing frequency. Then, it goes to quiescent state with a lower fractional exponent α = 0.75 , which is in good agreement with analytical results (see Fig. 3a-e). Next, with parameter set II, the integer order single neuron shows fast spiking while the input stimulus is on I = 45 . With the decrease of α = 0.84 , the firings transform into regular bursting, then with α = 0.82 and 0.80, it produces bursting however, the firing frequency decreases and more burst produces. Finally, it switches to quiescent state at α = 0.75 (see Fig. 3f-j).
Class II excitable neurons cannot generate low-frequency spikes. They are either in quiescent states or fire a train of spikes with larger frequency by a strong input current. The single M-L neuron with parameter set III shows fast spiking with α = 1 for I = 100 . With the decrease of α = 0.86 and 0.85, it generates MMBOs and MMOs. The firings switch to regular MMOs with further decrease of α = 0.84 , however it shows MMOs with lower firing frequency, i.e., the inter spike interval increases. Then, it goes to quiescent state α = 0.81 , i.e., converges to the fixed point of the system (see Fig. 3k-o). Now, we extend our study with the excitable slow-fast single 3D M-L neuron model (3) in the fractional domain with various parameter regimes that generate different bursting features, i.e., the number of spikes in each burst varies with diverse small and large amplitudes. With parameter set I, the single M-L model at α = 1 produces bursting with several number of spikes in each burst, however with the decrease of α = 0.9 and 0.8, the firing frequency decreases with longer time period, i.e., interspike interval increases between two burst and the amplitude of each spike decreases in the simultaneous bursting. For this parameter set, the unique fixed point is a saddle point of index two. It is observed that spike frequency adaptation occurs with the decrease of fractional-order exponents. Then it generates more spike frequency adaptation with further decrease of α = 0.7 (see Fig. 4a-d). Similarly, with parameter set II, the classical single neuron model shows bursting. With decreasing α = 0.95 and 0.75, it shows various bursting and then spiking behavior is observed with spike frequency adaptation and first spike latency at α = 0.5 (see Fig. 4e-h). Finally, for set III, the single neuron changes it behavior www.nature.com/scientificreports/ from bursting to fast spiking while α changes from one to α = 0.98 and 0.8. It switches to stable steady state with α = 0.6 , i.e., it converges to the locally asymptotically fixed point of the system (see Fig. 4i-l).

Network analysis of the fractional-order 2D M-L model
We investigate various dynamics of fractional-order M-L model in a random network architecture, where all the neurons are connected randomly with each other and with connection probability p . We construct an Erdös-Rényi network 56,57 of N = 100 M-L oscillators that is considered with mean node-degree �k� ∼ 7 for numerical simulations. The elaborate discussion of the network architecture is explained in the following subsections.

Dynamics of the network with two subpopulations.
To capture different firing activities of the network, first we consider a heterogeneous network of two different subpopulations depending on fractional-order (α) . Further, we consider that the fractional-order M-L neurons are electrically coupled through first state variable (u). The dynamics of network is studied using the following mathematical model The electrical coupling of the network is given by g e > 0 . The connection matrix of the network is represented by M = c ij N×N . Further, we divide the population of size N into two subpopulations depending on fractionalorder as m number of nodes have identical fractional-order α , reflecting oscillatory behavior and the remaining n nodes have fractional-order, β , reflecting excitable behavior. Thus, the total population size N can be expressed as N = m + n . First, we study the behavior of randomly connected class I excitable M-L neurons with two fractional-order exponents i.e., α = 1 and β = 0.75 . The total number of nodes in the network is N = 100 , and m = 60 & n = 40 i.e., we consider the network of N neurons with 60% oscillatory and 40% excitable neurons. In  www.nature.com/scientificreports/ the absence of coupling ( g e = 0 ), each oscillatory neuron in the network shows tonic spiking and the remaining neurons stay in quiescent state. With small increase in the electrical coupling g e = 0.0001 , oscillatory subpopulation still remains in tonic spiking mode and another subpopulation shows quiescent state. The time signals of two randomly connected nodes from two subpopulations are marked with red and blue lines (see Fig. 5a). The red signal is randomly chosen from the quiescent nodes and blue signal from the spiking nodes. The spatiotemporal plot reveals that the spiking nodes  are desynchronized to each other (Fig. 5e). The system behavior changes if we increase the coupling 100 folds ( g e = 0.01 ). Now, the subpopulation which was in quiescent state www.nature.com/scientificreports/ starts to exhibit bursting dynamics. Notably, the time interval of each burst is not periodic. Another subpopulation shows desynchronized irregular tonic spiking (see Fig. 5b, f). It is clear, if the coupling is increased in the mixed population, the periodic as well as quiescent nature vanishes and irregular bursting or spiking appears in the network. The entire network shows bursting dynamics with finite number of spikes in each burst with small increase of g e = 0.08 (see Fig. 5c, g). Here, two clusters with different amplitudes but the same phases are generated. Finally, at g e = 1 , the coupled network exhibits almost synchronized behavior by changing the firing activity to tonic spiking (Fig. 5d, h). Next, we increase the current stimulus I = 45 and the oscillatory subpopulation shows fast spiking and the other subpopulation remains in quiescent state. At weak coupling ( g e = 0.0001 ), the behavior of both subpopulations do not change. With the increase of coupling g e = 0.05 and 0.08, both the subpopulations start firing and show bursting dynamics. Finally, the activity of the random network changes to synchronized periodic spiking at g e = 1 (Fig. 5i-p). Next, we study the class II excitable M-L neurons in the same network architecture, i.e., one subpopulation is in quiescent state ( β = 0.81 ) and another subpopulation shows MMOs ( α = 0.86 ) in the absence of coupling. At weak coupling, the neurons are in desynchronized state ( g e = 0.0001, 0.01 ). Further increase of coupling, all the nodes in the network show spike frequency adaptation ( g e = 0.5 and g e = 1 ). It is clear that each of the two sets of oscillators exhibits almost complete synchronization showing similar type of bursting with identical phases and amplitudes (see Fig. 5q-x).
The reduced-order model with two subpopulations. In this section, a general and low-dimensional model description have been described to show that the reduced network exhibits the same feature as observed in the complete random network. In the intermediate coupling, we have observed a two cluster synchronization state that appears in the system. Motivated by this fact, we can write u 1 = u 2 = u 3 = · · · = u 60 = u α , and u 61 = u 62 = · · · = u 100 = u β . As we have considered the Erdös-Rényi graph, we can approximate the degree of each node/neuron by the average degree of the considered network [56][57][58][59][60][61][62] . Therefore, we can assume the degree of the j node k j = �k� . The number of spiking oscillators in the neighborhood of each oscillator is expected to be (1 − p e )k = p o k and the value will be p e k for the quiescent oscillators. Therefore, we can reconstruct a reducedorder model with two oscillators as follows www.nature.com/scientificreports/ where p e = n N and p o = m N represent the probabilities of excitable and oscillatory neurons respectively. We perform numerical simulations with the reduced order two coupled models, in which each subpopulation is captured by the identical fractional-order exponent exhibiting synchronous behavior for class I and class II M-L models. The numerical results suggest that the dynamics of the reduced-order model follows the same pattern as the entire graph when the two subpopulations follow cluster synchronization (see Fig. 6). For instance, bursting with two spikes can emerge for intermediate coupling in class I excitable system (Fig. 5k) separated by two clusters in the full network. Similar firing pattern exists in the reduced order model (Fig. 6g).

Conclusions
In this study, first we focused on the single neuron's dynamics that can change their response to various input statistics depending on the fractional exponents. We discussed the excitabilities of class I and class II conductance-based 2D M-L neurons that can be captured in realistic cortical neurons experimentally 2 . We examined corresponding bifurcation analysis of the classical-order model. The fractional exponent can trigger the firing variations that cannot be observed in the classical-order dynamics. The ability of the 2D M-L model for diverse spike responses including MMOs and MMBOs with fractional derivative is explored with spike frequency adaptation while decreasing the fractional exponent. The results demonstrate that the slow-fast 3D M-L in fractional domain provides various bursting patterns. It also changes its behavior from bursting to single train of spikes and fast spiking to irregular bursting for different suitable set of parameters that can not be captured in the classicalorder model for a fixed set of parameters. The FOD shows an alternative representation of the spiking-bursting responses with the changes of fractional exponents in the dynamical models. The approach summarized the multiple behavior of the single excitable model to certain stimulus variance with a memory dependent activities. We investigated the role of electrical coupling in a random network, where certain fraction of nodes are in quiescent states. If the oscillatory nodes are in MMBOs, the entire population would exhibit spike frequency adaptation. On the other hand, if the uncoupled oscillatory nodes are kept at fast tonic spiking zone, the entire population split into two clusters revealing periodic bursting in intermediate coupling and at higher coupling all the nodes show tonic spiking. Motivated by the cluster synchronization phenomena, we were also able to reduce the network into two coupled dynamics, which successfully captured the dynamics of the entire network during cluster synchronization.
We established distinct effects on different membrane voltage features considering the fractional-order derivative. We showed that the differences in the neuronal characteristics are due to the memory effects. Fractionalorder derivative provides rich dynamics and it is possible to explore realistic phenomena. These results demonstrate that the model and network provide a tractable approach to examine neuronal dynamics.

Data availibility
All numerically simulated data generated or analysed during this study are included in this submitted article.

Appendix: Numerical methods
In what follows, we present the numerical scheme used in our numerical simulations to approximate the solutions of the system of fractional-order differential equations of a variable X ≡ (u, v) T described by where α ∈ (0, 1) , f ≡ (h 1 , h 2 ) T and the Caputo-type derivative given by The Caputo-type fractional-order derivative is consistent with the physical initial and boundary conditions. In this case, the firing characteristics of the system are strongly independent of the initial conditions 4,12 . We can define, the initial conditions for the fractional-order system that can be handled using an analogy with the classical-order derivative. It includes a memory effect with a convolution between the classical-order derivative and a power of time. It is efficient to integrate all the previous activities of the function weighted by a function that follows power law dynamics. The numerical simulations of the systems are performed using the time step t = 0.1 . The simulation results with a smaller time step do not show any significant differences. The bifurcation diagrams of the fixed points of the dynamical model were computed using the MatCont software package 63 . We simulated the model (Eq. 3) using the L1 scheme 4,12,18 and approximated the fractional-order derivative as follows: (12) C d α uα dt α = −0.5g Ca (u α − V Ca )((1 + tanh(u α − V 1 ))/V 2 ) − v α g K (u α − V K ) − g L (u α − V L ) + I + g e p e (u β − u α ), (16) u(t N ) ≈ (dt) α Ŵ(2 − α)h 1 (u(t), v(t)) + u(t N−1 )